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Many  industrial  and  natural  processes,  like  nucleation  or  combustion  are  the  result  of  thousands 
of  interdependent  reactions.  Seemingly,  a  model  of  these  processes  demands  all  these  reactions 
to  be  taken  into  account  without  regard  to  their  ultimate  importance.  However,  despite  the  fast 
growth  of  computational  resources  of  the  last  few  decades,  this  extensive  modeling  remains  a 
daunting  task  except  when  massive  simplifications  are  made.  This  situation  mostly  originates 
from  two  problems:  high  computational  cost  of  the  methods  that  need  to  be  employed  to 
accurately  model  each  reaction  [1],  and  from  the  difficulty  of  finding  and  systematically 
exploring  all  the  possible  reactive  pathways.  Moreover,  even  if  issues  related  to  computational 
cost  will  be  resolved,  the  amount  of  information  needed  to  fully  describe  an  extensive  reaction 
network,  would  still  not  help  in  the  understanding  of  the  underlying  process  as  such  load  of 
information  would  be  difficult  to  interpret. 

Luckily,  in  many  cases  not  all  the  possible  reactions  have  the  same  importance,  as  often 
many  reaction  pathways  contribute  only  marginally  to  the  products’  formation.  This  fact  has 
been  leveraged  in  the  past  to  build  simplified  models,  often  called  reduced  mechanism  (RM), 
which  capture  the  key  aspects  of  a  more  detailed  description  while  making  it  more  tractable 
[2].  A  RM  is  built  on  two  distinct  components,  a  subset  of  reactions  and  their  corresponding 
rates;  it  is  generally  constructed  through  a  top-down  approach,  where  the  complete  reaction 
network  is  reduced  to  a  more  manageable  subset  [3].  However  except  for  the  most  simple 
cases,  key  reactions  are  difficult  to  identify  and  important  simplifications  are  made  in  the 
calculations  of  the  rates  [4],  which  contrasts  with  the  high-level  accuracy  required  to  correctly 
model  chemical  reactions  [5,  6].  Analogy  and  intuition  can  be  used  to  guess  which  are  the 
most  critical  reactions  ans  species  in  a  complex  reaction  network,  but  even  though  these  ideas 
are  then  verified  by  rigorous  testing,  this  approach  is  inefficient  as  it  relies  on  scientific  and 
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Chapter  1.  Introduction 


personal  bias. 

The  high  sensitivity  of  the  RM  to  the  methods  used  to  compute  the  rates  k  stems  from  the 
exponential  dependency  from  the  free  energy  (FE)  difference  between  products  and  transition 
state  A [7]: 


k  =  K-^—  exp(— jSAG*),  (1.1) 

h 

where  K  is  the  transmission  coefficient,  kg  is  the  Boltzmann  constant,  7  is  the  temperature, 
h  is  the  Planck’s  constant  and  /3  =  (kgT)~l.  From  Eq.  1.1  follows  that  any  error  in  computing 
AG  ‘:  is  likely  to  have  important  effects  on  the  values  of  the  derived  rates,  which  may  have 
catastrophic  repercussions  if  the  error  affects  an  early  branching  of  a  reaction  network. 

For  these  reasons,  instead  of  using  one  of  the  several  existing  methods  developed  to 
compute  the  reaction  dynamics  [8,  9,  10,  11,  12]  on  all  the  possible  reactions,  we  developed  a 
new  technique  to  identify  and  select  the  most  frequent  pathways  of  a  reaction  network  with 
minimal  computational  effort  by  using  an  acceleration-detection  scheme.  Briefly,  this  new 
approach  finds  the  critical  reactions  of  the  RM  gradually  by  iteratively  identifying  the  most 
common  products  of  given  reactants.  At  each  iteration,  the  products  of  the  previous  step  form 
the  pool  for  the  reactants  together  with  the  most  common  gas  species  present  under  specific 
conditions.  For  given  reactants,  the  most  frequent  reactions  are  obtained  by  accelerating  their 
dynamics  by  repeatedly  adding  small  amount  of  bias  to  the  FE  hyper-surface,  in  a  manner 
that  emulates  the  energy  transfer  associated  with  gas  phase  collisions.  As  a  result  a  list  of  the 
primary  pathways  is  produced;  the  corresponding  rates  can  then  to  be  computed  with  accurate 
ab-initio  techniques  in  order  to  build  the  final  RM. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


2. 1  Overview 

For  the  reasons  described  in  the  introductions,  instead  of  using  one  of  the  several  existing 
methods  developed  to  compute  the  reaction  dynamics  [8,  9,  10,  11,  12]  on  all  the  possible 
reactions,  we  developed  a  new  technique  to  identify  and  select  the  most  frequent  pathways  of  a 
reaction  network  with  minimal  computational  effort,  leaving  the  a  more  accurate  estimation  of 
the  rates  of  the  selected  reactions  to  high  level  ab-initio  calculations.  To  develop  our  approach, 
we  started  from  the  observation  that  for  any  given  reactant(s),  all  the  reaction  pathways  and 
rates  can  be  recovered  by  simply  observing  the  behavior  of  a  large  number  of  replicas  of 
the  same  system  for  a  long  time  and  counting  the  occurrences  of  each  reaction.  However, 
this  method  is  not  practical  because  it  requires  for  each  reactant  several  hundred  very  long 
simulations,  due  to  the  high  energy  barriers  commonly  involved  in  chemical  reactions.  To 
make  this  idea  applicable,  we  employed  an  acceleration-detection  scheme,  where  the  dynamics 
of  all  the  system’s  replicas  are  accelerated  until  a  reaction  is  detected.  The  simulations  are  then 
interrupted  and,  after  all  the  replicas  are  terminated,  the  frequency  of  each  reactive  pathway  is 
calculated.  This  method  has  several  advantages: 

•  it  does  not  require  a  priori  knowledge  of  the  pathways  or  transition  states; 

•  it  does  not  rely  on  the  life  time  or  stability  of  the  products,  as  the  simulations  are 
interrupted  as  soon  as  the  reaction  happens,  and  therefore  this  approach  can  also  handle 
pathways  where  chemical  activation  plays  an  important  role. 

For  given  reactant(s)  this  approach  will  produce  a  list  of  pathways  and  by  repeating  the 
same  procedure  only  for  the  most  frequently  observed  reactions,  the  key  pathways  of  the  entire 
reaction  network  are  obtained  without  the  need  to  either  map  the  complete  reaction  network  or 
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Chapter  2.  Methodology 


arbitrarily  select  pathways. 


2.2  Acceleration 

Since  the  idea  of  accelerating  systems  to  overcome  energetic  barriers  is  not  new  [13],  we 
used  the  basic  algorithm  that  is  use  in  Metadynamics  (META)  [14],  an  already  well-tested 
method  [15]  that  uses  a  history-dependent  bias  to  favor  the  exploration  of  new  states.  Briefly, 
the  well-tempered  Metadynamics  (WtM)  technique  was  introduced  to  reconstruct  the  FE 
landscape  of  a  system  projected  on  a  few  coordinates,  called  collective  variables  (CVs),  by 
adding  an  history  dependent  bias  potential  that  favors  the  exploration  of  configuration  that  are 
far  from  equilibrium.  Moreover,  to  get  the  added  potential  to  convergence  to  the  underlying 
FE,  in  the  WtM  algorithm  the  amount  of  added  bias  is  gradually  reduced  in  previously  visited 
values  of  the  CVs.  The  amount  of  this  reduction  is  controlled  by  the  bias  factor  (bf),  y 


7 — 


T+AT 

T 


(2.1) 


where  T  is  the  temperature  of  the  system  and  T  +  AT  is  the  fictitious  temperature  of  the 
CVs,  as  in  the  long  time  limit  the  probability  distribution  of  the  CVs  is 


FE 

P( CV)  oc  e  vkbT  (2.2) 

The  reader  can  find  a  more  complete  description  of  META  and  wtM  in  the  cited  papers. 

Even  though  META  was  originally  introduced  to  reconstruct  FE  landscapes,  here  we  employ 
it  only  to  accelerate  the  system  dynamics,  simplifying  the  requirements  on  the  definition  of 
the  CV.  In  particular,  we  use  the  algorithm  to  bias  the  potential  energy,  forcing  the  reacting 
molecule  to  experience  energy  fluctuations  typical  of  higher  temperatures  [16].  This  approach 
can  be  viewed  as  a  way  to  force  the  molecule  under  investigation  to  experience  several 
collisions  with  virtual  particles  that  only  increase  its  internal  energy,  effectively  accelerating 
its  reactions  in  the  high-pressure  regime. 

With  the  addition  of  the  potential  energy  bias  <P,  the  effective  FE  experienced  by  the 
system  becomes  AG)  +  <T>,  which  substituted  in  Eq.  1.1  can  be  used  to  compute  the  probability 
Pi  to  observe  a  specific  reaction  i: 


KjMi  exp  ( — /3  AG  * ) 

Pi  —  - -r-  (2.3) 

ZjKjMjcxp(-PAGj) 

where  M,  is  the  pathway  degeneracy  of  the  z'-th  reaction,  and  the  summation  is  performed 
over  all  the  possible  reactions. 

While  biasing  the  potential  energy  is  an  effective  and  general  way  to  accelerate  the  system 
reactivity  [17],  at  the  same  time  it  is  not  an  appropriate  quantity  to  use  to  monitor  the  evolution 
of  the  reactions,  since  in  many  cases  it  is  unable  to  distinguish  between  products  and  reactants. 
A  more  apt  choice  for  the  reaction  detection  is  the  measure  of  the  molecular  connectivity,  like 
the  recently  introduced  Social  PeRmutation  INvarianT  coordinates  (SPRINT)  [18].  This  class 
of  reaction  coordinates  defines  the  connectivity  of  a  system  of  N  atoms  with  a  N-dimensional 
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2.3  SPRINT 
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vector  by  using  spectral  graph  theory  and  including  both  local  and  long-range  system  topology 
information.  SPRINT  coordinates  have  already  been  successfully  used  to  differentiate  and 
cluster  molecular  structures  [19],  by  considering  the  evolution  of  the  Euclidean  norm  of  the 
difference  between  instantaneous  and  average  value  of  the  SPRINT. 


2.3  SPRINT 

The  SPRINT  CVs  are  related  to  the  concept  of  coordination  number  (CN)  of  each  atom.  As  the 
CN  is  so  central  for  the  detection  part  of  the  algorithm  to  work  it  is  worth  tune  its  parameters  so 
that  they  are  optimized  for  the  Fast  Exploration  of  Reaction  Network  (FERN)  algorithm.  Due 
to  the  derivability  requirements  off  all  the  CVs  in  the  META  technique,  a  smooth  asymmetric 
step  function  is  needed.  We  employed  the  one  described  by  the  following  equation  (m  >  n ): 


cn(r ) 


r  <  d0 


r  >  do 


(2.4) 


To  define  this  function  four  parameters  needs  to  be  provided:  do,  fo,  n,  and  m.  While 
do  simply  shifts  the  curve  on  the  r  axis,  and  represents  the  maximum  distance  for  which  a 
specific  CN  is  always  equal  to  one,  the  role  of  the  other  parameters  is  slightly  more  complex. 
To  understand  their  effect  we  can  consider  the  three  scenarios: 


r  — >  do 

l-r 

(2.5) 

cn  ps - 

1  -  rm 

r  — >  inf 

cn  ps  r"-w 

(2.6) 

do  +  fo 

cn  —  nj  m 

(2.7) 

From  equation  2.5  we  can  see  that  the  behavior  a  short  distance  is  dictated  by  m,  while 
from  equation  2.6  can  be  deduced  that  in  order  to  have  a  slow  decaying  function  we  need 
to  keep  m  —  n  small;  finally  the  third  case  (equation  2.7)  shows  that  ro  measures  the  interval 
required  to  go  from  1  to  n/m. 

The  exact  values  to  assign  to  these  variables  depends  on  the  goal  of  the  simulations.  As 
we  are  interested  in  bond  breaking,  the  requirements  that  we  want  to  meet  are: 

•  minimize  the  effect  on  SPRINT  of  thermal  oscillation  for  a  given  configuration. 

•  avoid  the  collapse  of  SPRINT  values  to  0  that  follows  excessive  separation  of  atoms.  Is 
worth  nothing  that  this  requirement  is  only  relatively  important  since  each  reaction  is 
interrupted  as  soon  as  a  clear  change  in  the  connectivity  is  observed  and  therefore,  a 
weak  connectivity  even  after  a  bond  is  broken  is  not  strictly  required. 

•  assign  the  same  (integer)  value  of  n  and  m  among  all  the  all  the  types  of  bonds,  due 
to  programming  restraints.  While  this  requirement  can  be  remove  by  rewriting  the 
implementation  of  the  SPRINT  calculation,  it  was  not  removed  as  it  is  of  no  influence  on 
the  results. 
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With  these  requirements  in  mind,  we  first  chose  n  —  2  and  m  =  3,  which  give  the  least 
steep  function  with  the  slowest  decay  (n  —  1  causes  a  discontinuity  in  the  first  derivative  at 
do).  Then  we  fine  tuned  the  values  of  do  and  ro  for  each  bond  type  by  computing  the  FE  as  a 
function  of  the  bond  length.  As  our  test  are  performed  on  hydrocarbons  (see  next  chapter), 
here  we  report  the  results  for  the  C— C,  C— H,  and  H— H  bonds,  but  the  procedure  is  equivalent 
for  other  type  of  atoms. 

2.3.1  C-C,  C-H,  and  H-H  bond  analysis 

The  FE  were  computed  by  using  wtM  biased  molecular  dynamics  (md)  simulations  performed 
with  the  LAMMPS  [20]  software  coupled  with  the  PLUMED  [21]  plugin  (version  1.3).  The 
reaction  were  carried  by  employing  the  AIREBO  [22]  classic  reactive  force  field  (FF).  The 
equations  of  motion  were  integrated  with  a  timestep  of  0.1  fs;  systems  were  simulated  for 
50  ns,  at  a  temperature  of  2500  K  maintained  with  a  Langevin  thermostat  [23]  with  a  time 
constant  X[an  =  5  fs;  the  cutoff  distance  multiplier  was  set  to  3.  As  the  Gaussian  shaped  bias 
with  an  initial  height  h  of  0.02  eV  and  a  width  a  of  0.025  nm  was  deposited  every  50  fs  on  the 
specific  bond  distance  while  employing  a  BF  of  6.  The  bond  exploration  of  the  phase  space 
was  controlled  by  placing  an  harmonic  "wall"  at  a  distance  of  0.5  nm  with  an  elastic  constant 
of  50  eV/nm2 
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Figure  2.1:  (right  axis)  FE  as  a  function  of  the  H— H  bond  length  in  H2.  (left  axis)  CN  as  a 
function  of  H— H  distance.  The  gray  vertical  line  indicates  the  approximate  location  of  the 
transition  state. 


Due  to  the  cross-dependence  of  the  parameters  for  the  different  bonds  on  the  value  of  the 
SPRINT  we  first  defined  the  CN  parameters  for  H— H  bond  by  studying  H2  (Figure  2.1),  we 
then  moved  to  the  C— H  bond  by  looking  into  CH4  (Figure  2.2),  and  finally  on  the  C— C  bond 
by  analyzing  n-C7H16  (Figure  2.2). 
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Figure  2.2:  (right  axis)  FE  as  a  function  of  the  C— H  bond  length  in  CH4.  (left  axis)  CN  as  a 
function  of  C— H  distance.  The  gray  vertical  line  indicates  the  approximate  location  of  the 
transition  state. 
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Figure  2.3:  (right  axis)  FE  as  a  function  of  the  C— C  bond  length  of  two  terminal  carbons  in 
n-heptane.  (left  axis)  CN  as  a  function  of  C— C  distance.  The  gray  vertical  line  indicates  the 
approximate  location  of  the  transition  state. 


By  analyzing  the  FEs,  we  found  that  by  assigning  ro  =  0.05  and  a  value  to  do  0.03  nm 
smaller  than  the  position  of  the  transition  state  (shown  with  a  gray  vertical  line  in  the  figure), 
allows  to  detect  only  configuration  which  have  enough  energy  to  potentially  cross  the  FE 
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barrier,  while  at  the  same  time  ignoring  "unproductive"  thermal  fluctuations  even  at  high 
temperature.  The  final  values  of  do  for  the  H— H,  C— H,  and  C— C  bonds  are  0.14  nm,  0.15 
nm,  and  0.182  nm,  respectively. 
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To  illustrate  the  salient  details  of  our  acceleration-detection  approach,  we  applied  our  method¬ 
ology  to  three  different  systems  of  increasing  complexity.  First,  we  studied  the  behavior  of  a 
monatomic  particle  subjected  to  an  analytical  two  dimensional  FE  potential.  This  is  a  simple 
system  for  which  all  the  properties  are  known  exactly  and  therefore  can  be  used  to  test  the 
correctness  of  the  method’s  foundations.  The  second  case  we  considered  is  the  first  step  of 
the  unimolecular  dissociation  of  ethane.  Compared  to  the  first,  this  system  has  the  added 
complexity  of  the  rotational  and  vibrational  modes  of  the  molecules,  so  we  could  test  and 
discuss  all  the  issues  related  to  the  thermostatting  and  temperature  control.  Moreover,  as  we 
are  employing  a  polyatomic  systems,  we  illustrate  the  how  to  set  up  the  parameters  used  to 
compute  the  SPRINT.  The  third  and  last  system,  we  studied  is  the  reactivity  of  /-decalin  with 
methyl  radical.  Not  only  the  number  of  pathways  for  this  system  is  larger,  but  also  details  on 
bimolecular  reactions  are  presented. 


3. 1  Monoatomic  particle  in  2D  potential 

The  behavior  of  a  single  particle  in  a  customizable  analytic  potential  is  an  excellent  test  for 
the  performance  of  the  FERN  method.  The  potential  was  shaped  to  reproduce  three  minima 
separated  by  two  barriers  (symmetric  or  asymmetric,  depending  on  the  values  of  the  chosen 
parameters),  as  shown  in  Figure  3.1. 

The  potentials  are  described  by  the  following  equation: 

U(x,y )  =  axx 6  +  7  axx4  +  1 2  axx2  +  dxx  +  ayy2  (3.1) 

where  dx  can  be  used  to  change  the  difference  in  energy  between  the  barriers.  All  the 
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Figure  3.1:  Example  of  the  two  dimensional  potential  experienced  by  the  single  atom  ( ax  —  3, 
dx  —  3.3,  ay  —  3).  The  energy  is  expressed  in  kcal/mol  and  the  contour  lines  are  separated  by 
5  energy  units. 


simulations  were  started  from  the  central  basin  and  the  rates  of  the  formation  of  “products” 
were  computed  by  monitoring  the  position  of  the  particles  on  the  x  axis,  for  in  this  simple  case 
there  is  no  chemical  connectivity  that  justifies  the  use  of  SPRINT.  The  simulations  were  as 
soon  as  the  atom  entered  either  the  negative  (x  <  —1)  or  the  positive  (x  >  1)  basins. 

By  tweaking  the  parameters  that  define  the  potential,  we  created  several  scenarios  with 
different  barrier  heights  (from  about  4  to  80  times  kgT)  and  different  degrees  of  asymmetry 
between  the  negative  and  positive  basins  (with  a  difference  between  the  two  barriers  ranging 
from  0  to  about  15  kgT).  We  analyzed  the  effect  of  the  thermostat  and  a  variety  of  META 
parameters,  in  particular  deposition  rate  ( 8h ),  bias  shape  (height  hn  and  width  Oh  of  the 
Gaussian  shaped  bias)  and  BF,  when  employing  the  well-tempered  version  [11]. 

Initial  tests  were  performed  on  symmetric  systems  ( dx  —  0)  with  (Table  3.2)  and  without 
(Table  3.1)  a  thermostat.  Both  wtM,  with  different  values  of  the  BF  and  "plain"  META  (BF 
=  co)  were  employed.  For  each  system  we  reported  the  average  and  standard  deviation  (over 
all  the  runs)  of  the  temperature)  (7)  and  Oj ,  both  in  K),  and  of  the  simulation  time  ((/)  and  ot, 
both  in  ps). 

In  these  tests  we  expect  to  statistically  observe  the  atom  end  in  the  positive  and  negative 
well  the  same  number  of  time,  therefore  these  tests  can  be  used  to  evaluate  the  dependency  of 
the  method  on  different  simulations  parameters  that  have  nothing  to  do  with  the  underlying 
physical  system.  In  all  cases,  the  discrepancy  between  the  observed  and  the  expected  probabil¬ 
ity  are  within  a  95%  confidence  interval.  Interestingly  however,  increasing  the  amount  of  bias 
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Table  3.1:  Results  for  different  symmetric  (dx  —  0)  two  dimensional  potentials  systems 
simulated  without  a  thermostat.  Bias  shape  parameters  (hfj ,  Oh)  as  well  as  FE  barriers  (AGCC, 
A G^ch)  are  expressed  in  kcal/mol,  while  bias  addition  frequency  <5#  is  expressed  in  ps;  masses 
are  listed  in  Dalton.  Temperature  average  and  standard  deviation  ((T),  Of)  are  expressed  in  K; 
time  average  and  standard  deviation  (( t ),  ot)  are  in  ps.  The  "Neg."  and  "Pos."  columns  refer 
to  the  number  of  observed  simulation  ending  in  the  negative  (x  <  —1)  or  the  positive  (x  >  1 ) 
basins,  respectively. 


BF 

hH 

oH 

8h 

&X 

dy 

Runs 

(T) 

Of 

(0 

Ot 

Neg. 

Pos. 

OO 

0.05 

0.05 

1 

2 

2 

60 

639 

300 

2.19 

2.05 

31 

29 

3 

0.05 

0.05 

1 

2 

2 

60 

570 

246 

2.77 

2.89 

32 

28 

3 

0.01 

0.01 

1 

2 

2 

60 

1017 

843 

2.08 

1.87 

28 

32 

3 

0.01 

0.01 

0.1 

2 

2 

60 

1004 

288 

1.01 

1.31 

35 

25 

3 

0.01 

0.01 

0.05 

2 

2 

60 

965 

336 

0.74 

0.69 

29 

31 

3 

0.1 

0.1 

10 

2 

2 

60 

301 

112 

24.10 

22.55 

25 

35 

3 

0.01 

0.03 

0.05 

2 

2 

60 

867 

255 

1.58 

1.38 

34 

26 

3 

0.01 

0.1 

0.05 

2 

2 

60 

403 

155 

10.63 

4.82 

26 

34 

10 

0.01 

0.1 

0.05 

3 

3 

100 

1250 

420 

2.17 

2.34 

47 

53 

3 

0.05 

0.05 

1 

4 

2 

60 

1286 

466 

2.39 

2.78 

32 

28 

3 

0.05 

0.05 

1 

4 

4 

60 

1341 

541 

1.91 

1.93 

27 

33 

10 

0.01 

0.1 

0.05 

10 

2 

100 

3587 

1237 

2.09 

2.61 

58 

42 

3 

0.01 

0.1 

0.05 

10 

10 

160 

3507 

1113 

2.36 

1.90 

83 

77 

10 

0.01 

0.1 

0.05 

10 

10 

160 

3823 

1372 

2.70 

2.28 

79 

81 

25 

0.01 

0.1 

0.05 

10 

10 

160 

3980 

1264 

1.65 

2.26 

89 

71 

50 

0.01 

0.1 

0.05 

10 

10 

160 

3724 

1203 

1.90 

2.28 

90 

70 

100 

0.01 

0.1 

0.05 

10 

10 

160 

3904 

1392 

1.80 

2.12 

86 

74 

deposited  at  each  interval  does  not  lead  to  shorter  simulations.  As  expected  in  the  simulations 
with  the  thermostat,  we  found  that  very  small  T  are  almost  never  beneficial,  slowing  the 
diffusion  of  the  atom  on  the  FE  surface. 

A  second  set  of  tests  was  performed  with  asymmetric  potentials,  again  with  (Table  3.3)  and 
without  (Table  3.4)  the  presence  of  a  thermostat.  As  in  the  previous  case,  the  application  of  a 
thermostat  is  expected  not  to  affect  the  ratio  between  the  rates  due  to  the  lack  of  vibrational 
degrees  of  freedom  of  the  system.  The  AG':  needed  to  compute  the  theoretical  ratio  was 
evaluated  directly  from  3.1,  since  there  is  no  relevant  difference  in  the  entropic  contribution  of 
the  two  pathways. 

A  subset  of  the  results  (selected  for  clarity)  is  shown  in  Figure  3.2. 

In  all  cases  we  found  an  excellent  agreement  between  the  pathway  probability  computed 
with  our  approach  and  the  one  predicted  by  using  the  analytical  values  of  AG  :  in  Equation  2.3. 
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Table  3.2:  Results  for  different  symmetric  two  dimensional  potential  (ax  =  4,  ay  =  2,  dx  —  0) 
systems  simulated  with  a  thermostat.  Bias  shape  parameters  (hn,  Oh)  as  well  as  FE  barriers 
(AGf-c,  A G^ch)  are  expressed  in  kcal/mol,  while  bias  addition  frequency  8h  is  expressed  in 
ps;  masses  are  listed  in  Dalton.  Temperature  average  and  standard  deviation  ((T),  Oj)  are 
expressed  in  K;  time  average  and  standard  deviation  ((f),  o,)  are  in  ps.  The  "Thermo"  column 
refers  to  the  type  of  thermostat  (Lan  for  Langevin  and  NHclO  for  a  Nose-Hoover  chain  of 
length  10),  applied  with  a  time  constant  r  (in  ps).  The  "Neg."  and  "Pos."  columns  refer  to  the 
number  of  observed  simulation  ending  in  the  negative  {x  <  —  1)  or  the  positive  (x  >  1)  basins, 
respectively. 


BF 

hH 

oH 

8h 

Thermo 

T 

Runs 

(T) 

Ot 

(0 

ot 

Neg. 

Pos. 

3 

0.05 

0.05 

1 

Lan 

10 

60 

1462 

427 

1.95 

2.72 

35 

25 

3 

0.05 

0.05 

1 

Lan 

1 

60 

1575 

430 

1.04 

0.72 

34 

26 

3 

0.05 

0.05 

1 

Lan 

0.1 

58 

1944 

644 

3.86 

9.48 

30 

28 

3 

0.05 

0.05 

1 

NHclO 

10 

60 

1518 

425 

1.02 

2.04 

26 

34 

3 

0.05 

0.05 

1 

NHclO 

1 

59 

1535 

550 

1.57 

6.44 

26 

33 

3 

0.05 

0.05 

1 

NHclO 

0.1 

60 

1453 

451 

7.22 

12.56 

31 

29 

3 

0.01 

0.1 

0.05 

NHclO 

0.1 

38 

901 

128 

21.30 

11.61 

21 

17 

3 

0.01 

0.1 

0.05 

NHclO 

10 

100 

897 

440 

4.66 

3.82 

45 

55 

Table  3.3:  Results  for  different  asymmetric  two  dimensional  potentials  systems  simulated 
with  a  thermostat.  The  FE  barriers  (AG)C,  AG^H)  are  expressed  in  kcal/mol,  while  masses 
are  listed  in  Dalton.  The  bias  (hn  =  0.01  kcal/mol,  Oh  =  0. 1  kcal/mol)  was  added  with  a 
frequency  8h  of  0.05  ps.  Temperature  average  and  standard  deviation  ( (T).  Or)  are  expressed 
in  K;  time  average  and  standard  deviation  ((f),  ot)  are  in  ps.  T  (in  ps)  is  the  time  constant 
of  the  Nose-Hoover  chain  thermostat  (length  10)  "Neg."  and  "Pos."  columns  refer  to  the 
number  of  observed  simulation  ending  in  the  negative  (x  <  —  1)  or  the  positive  (x  >  1)  basins, 
respectively. 


BF 

ax 

Cly 

d-x 

T 

Runs 

!T) 

Oj 

(0 

ot 

Neg. 

Pos. 

3 

2 

2 

0.6 

0.01 

120 

508 

19 

19.87 

14.24 

107 

13 

3 

2 

2 

0.6 

0.05 

52 

567 

61 

7.27 

3.77 

46 

6 

3 

2 

2 

0.6 

0.1 

98 

656 

153 

7.06 

10.86 

86 

12 

3 

2 

2 

0.6 

0.05 

100 

1389 

432 

1.17 

2.86 

64 

36 

1.5 

2 

2 

0.6 

0.05 

97 

1228 

275 

3.52 

3.91 

74 

23 

50 

2 

2 

0.6 

0.05 

100 

763 

215 

3.29 

5.67 

83 

17 

50 

8 

8 

0.15 

0.05 

65 

1657 

204 

17.95 

5.09 

47 

18 

50 

12 

12 

0.1 

0.05 

32 

2580 

771 

19.90 

15.09 

21 

11 
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Table  3.4:  Results  for  different  asymmetric  two  dimensional  potentials  systems  simulated 
without  a  thermostat.  Bias  shape  parameters  (hn,  Oh)  as  well  as  FE  barriers  (A GJ^C,  A GVCH) 
are  expressed  in  kcal/mol,  while  bias  addition  frequency  Sh  is  expressed  in  ps;  masses  are 
expressed  in  Dalton.  Temperature  average  and  standard  deviation  ({T),  oy)  are  expressed  in 
K;  time  average  and  standard  deviation  ((?),  or)  are  in  ps.  The  "Neg."  and  "Pos."  columns 
refer  to  the  number  of  observed  simulation  ending  in  the  negative  (x  <  —1)  or  the  positive 
(x  >  1)  basins,  respectively. 


BF 

hH 

oh 

SH 

ax 

ay 

d-x 

Runs 

(T) 

Oj 

<*> 

ot 

Neg. 

Pos. 

3 

0.01 

0.1 

0.05 

2 

2 

1 

60 

508 

193 

5 

5 

59 

1 

3 

0.05 

0.05 

1 

4 

4 

1 

60 

1327 

392 

2 

2 

56 

4 

3 

0.05 

0.05 

1 

3 

3 

1.4 

60 

940 

361 

2 

2 

56 

4 

CO 

0.05 

0.05 

1 

3 

3 

1.4 

60 

989 

495 

2 

3 

57 

3 

3 

0.05 

0.05 

1 

3 

3 

1.4 

60 

488 

286 

5 

4 

58 

2 

3 

0.01 

0.01 

0.1 

3 

3 

1.4 

60 

1164 

424 

1 

2 

56 

4 

3 

0.05 

0.05 

1 

3 

3 

1.4 

60 

1252 

297 

1 

2 

59 

1 

3 

0.01 

0.1 

0.05 

3 

3 

0.6 

60 

587 

328 

6 

3 

55 

5 

25 

0.01 

0.1 

0.05 

2 

2 

0.6 

97 

1227 

603 

4 

8 

76 

21 

3 

0.01 

0.1 

0.05 

2 

2 

0.6 

99 

747 

245 

3 

3 

88 

11 

3 

0.01 

0.1 

0.05 

10 

10 

0.6 

86 

1440 

975 

7 

4 

85 

1 

Figure  3.2:  Comparison  between  theoretical  (full  symbols)  and  computed  (empty  symbols) 
probability  for  the  transition  from  the  central  (—1  <  x  <  l)to  the  negative  basin  (x  <  —  1) 
for  different  2d  potentials.  Data  from  simulations  performed  with  (triangles)  and  without  a 
thermostat  (circles);  vertical  bars  show  the  95%  confidence  interval. 
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3.2  Ethane  decomposition 

As  a  second  example  we  considered  the  first  step  of  unimolecular  decomposition  of  ethane  at 
high  temperatures  (approximately  between  1000  and  3000  K).  As  with  the  previous  system 
only  two  pathways  are  possible,  a  C— C  or  a  C— H  bond  breaking,  but  this  time  the  reactions 
have  different  multiplicities  (one  and  six,  respectively)  and  dissimilar  entropic  contributions. 
For  all  simulations  we  used  adaptive  intermolecular  reactive  bond-ordered  potential  AIREBO 
[22].  While  classical  reactive  ff  are  not  necessarily  accurate  compared  to  ab-initio  or  density 
functional  theory  methods,  they  provides  a  consistent  and  computationally  light  framework 
to  test  our  method.  Since  the  method  itself  is  not  dependent  in  any  way  on  the  underlying 
potential,  this  choice  affects  the  result  but  not  the  validity  of  the  tests. 

3.2.1  fe  calculations 

In  order  to  validate  the  rates  computed  with  FERN  we  needed  an  estimate  of  the  rate  based  on 
the  FE  surface  at  different  temperatures  for  both  the  possible  reactions.  The  knowledge  of  the 
FE  difference  allows  to  use  Equation  2.3  to  compute  the  theoretical  rates  under  each  condition. 
To  obtain  the  values  of  A Gf  at  different  temperatures,  we  interpolated  the  FE  computed  at  500, 
1000,  1500  and  2000  K  and  assumed  the  transmission  rates  equal  for  all  the  reactions  (see 
Figures  3.3  and  3.4). 


Temperature  (K) 


Figure  3.3:  Computed  FE  for  the  C— C  bond  breaking  in  the  ethane  molecule  at  different 
temperatures.  Results  are  shown  with  the  95%  confidence  interval.  The  black  line  shows  the 
exponential  fitting. 
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Temperature  (K) 

Figure  3.4:  Computed  FE  for  the  C— H  bond  breaking  in  the  ethane  molecule  at  different 
temperatures.  Results  are  shown  with  the  95%  confidence  interval.  The  black  line  shows  the 
exponential  fitting. 


The  values  at  the  four  reference  temperature  were  computed  by  employing  WTM  and  the 
relevant  bond  distance  as  a  CV.  For  these  simulations  we  employed  MD  settings  similar  to  the 
one  reported  in  subsection  2.3.1,  with  the  only  exception  of  the  simulation  length  (150  ns)  and 
the  value  of  the  BF  which  was  modified  so  that  T *BF  was  almost  constant  (about  3500  K). 

The  data  points  reported  in  Figures  3.3  and  3.4  is  the  FE  difference  between  the  bonded 
and  the  transition  state.  The  first  is  defined  as  the  region  around  the  minimum  (MIN)  that  has 
up  to  2 kgT  more  energy  than  MIN,  while  the  latter  is  the  region  around  the  maximum  that  has 
up  to  \ksT  less  energy  than  the  barrier. 

3.2.2  Results 

The  FERN  method  was  tested  on  a  variety  of  different  systems,  related  to  both  the  simulation 
settings  and  the  physic  of  the  system  (isotopic  effect)  as  listed  in  Table  3.5 

We  did  not  apply  a  thermostat  as  its  efficiency  is  not  constant  for  all  the  frequencies 
and  therefore  its  use  can  radically  influences  the  results.  Instead,  the  control  on  the  final 
temperature  was  obtained  by  changing  the  biasing  parameters,  while  still  maintaining  each 
added  bias  relatively  small  (normally  less  than  1/1 0th  of  kgT  for  both  Gaussian  height  and 
width).  Figure  3.5  shows  a  selection  (chose  for  clarity)  of  the  results,  while  the  complete  list 
is  reported  in  Table  3.6. 

As  can  be  seen  in  Figure  3.5  the  predicted  probability  of  a  pathway  is  recovered  in  a  wide 
range  of  temperatures  and  simulations  parameters,  with  a  minimal  computational  cost.  Even 
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Table  3.5:  List  of  the  different  systems  used  to  test  the  ethane  reactivity.  Bias  shape  parameters 
(hfj ,  o h)  as  well  as  FE  barriers  (A G^c,  A G^CH)  are  expressed  in  kcal/mol,  while  bias  addition 
frequency  8h  is  expressed  in  ps;  masses  are  in  Dalton. 


Label 

BF 

hH 

SH 

mc 

mH 

^cc 

^g'ch 

Runs 

EtHl 

3 

0.01 

0.1 

0.05 

12 

1 

129 

119 

100 

EtH2 

OO 

0.01 

0.1 

0.05 

12 

1 

129 

118 

200 

EtH3 

3 

0.1 

0.25 

1 

12 

1 

136 

122 

39 

EtH4 

3 

0.25 

0.25 

1 

12 

1 

129 

120 

40 

EtH5 

3 

0.25 

0.25 

1 

12 

1 

132 

121 

40 

EtH6 

3 

1 

0.25 

1 

12 

1 

129 

120 

40 

EtH7 

3 

1 

0.25 

1 

12 

1 

129 

120 

40 

EtH8 

3 

1 

0.25 

1 

12 

1 

129 

120 

38 

EtH9 

3 

1 

0.25 

1 

12 

1 

129 

120 

39 

EtHIO 

3 

1 

0.25 

1 

12 

1 

132 

121 

40 

EtHll 

3 

1 

0.25 

1 

12 

1 

132 

121 

34 

EtHl  2 

250 

1 

0.25 

1 

12 

1 

129 

120 

40 

EtHl  3 

250 

1 

0.25 

1 

12 

1 

129 

120 

40 

EtHl  4 

250 

1 

0.25 

1 

12 

1 

129 

120 

40 

EtHl  5 

3 

4 

0.25 

1 

12 

1 

133 

121 

35 

EtHl  6 

OO 

0.01 

0.1 

0.05 

12 

2 

129 

119 

197 

EtHl  7 

OO 

0.01 

0.1 

0.05 

12 

3 

129 

119 

197 

EtHl  8 

3 

0.01 

0.1 

0.05 

12 

4 

131 

121 

70 

EtHl  9 

OO 

0.01 

0.1 

0.05 

12 

4 

129 

119 

194 

EtH20 

3 

1 

0.25 

1 

12 

4 

129 

120 

40 

EtH21 

3 

1 

0.25 

1 

12 

4 

147 

124 

29 

EtH22 

OO 

0.01 

0.1 

0.05 

12 

5 

130 

120 

182 

EtH23 

OO 

0.01 

0.1 

0.05 

12 

9 

131 

121 

91 

EtH24 

3 

1 

0.25 

1 

12 

9 

148 

124 

11 

without  any  specific  optimization,  the  average  simulation  time  is  most  of  the  time  10  ps  or 
less.  Considering  a  few  hundred  simulations  for  each  system,  the  total  required  simulation 
time  is  on  the  order  of  a  few  nanoseconds  and  can  be  further  tuned  by  modifying  the  number 
of  simulations  performed  on  each  node. 

Another  interesting  result  is  the  conservation  of  some  aspect  of  the  dynamics.  Despite  not 
being  in  any  part  of  the  FERN  formulation,  this  feature  of  can  be  observed  when  considering 
the  isotopic  effect  associated  with  the  change  in  mass  of  the  H.  As  AG^  is  not  affected  by 
the  change  of  mass,  the  number  of  observed  pathway  is  only  minimally  affected  with  respect 
to  the  inn  =  1  case,  but  the  isotopic  effect  is  recovered  in  the  average  reaction  time  {/),  which 
is  dominated  by  the  average  time  for  the  C-H  reaction.  Not  only  the  time  increases  with  the 
mass  but  also  the  correct  mass  dependence  (°c  yjtrin)  is  recovered  as  shown  in  Figure  3.6. 
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Figure  3.5:  ( lefty-axis )  Comparison  between  theoretical  (full  circles)  and  computed  (empty 
circles)  probability  for  the  C— H  bond  breaking.  Vertical  bars  show  the  95%  confidence 
interval,  {right  y-axis )  Squares  indicate  the  average  simulation  time  for  each  system;  each 
point  represents  the  average  of  40  to  200  simulations  (see  Table  3.5). 


Figure  3.6:  Average  reaction  time  for  the  C— H  bond  breaking  for  different  H  masses.  Data 
with  95  %  confidence  interval  are  shown  in  red;  fitting  function  plotted  in  black. 
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Table  3.6:  Results  for  the  simulations  of  the  ethane  reactivity.  Temperature  average  and 
standard  deviation  ((T),  aT)  are  expressed  in  K.  Time  average  and  standard  deviation  ((/),  ot) 
are  in  ps.  The  "CC"  and  "CH"  columns  refer  to  the  number  of  observed  simulation  ending 
with  a  C— C  or  a  C— H  bond  breaking. 


Label 

(T) 

Oj 

(0 

ot 

CC 

CH 

EtHl 

2460 

505 

7.79 

2.14 

2 

98 

EtH2 

2887 

654 

5.11 

2.93 

5 

195 

EtH3 

1057 

374 

47.10 

31.65 

1 

38 

EtH4 

2905 

697 

3.89 

2.16 

1 

39 

EtH5 

1417 

593 

13.05 

9.41 

0 

40 

EtH6 

4134 

970 

0.94 

0.51 

0 

40 

EtH7 

2738 

989 

1.26 

0.92 

2 

38 

EtH8 

2632 

826 

1.25 

0.97 

0 

38 

EtH9 

2760 

1021 

1.38 

0.94 

0 

39 

EtHIO 

1358 

593 

2.97 

1.72 

1 

39 

EtHll 

1336 

720 

2.96 

1.98 

0 

34 

EtHl  2 

3985 

799 

0.94 

0.49 

0 

40 

EtHl  3 

2634 

761 

1.45 

0.98 

1 

39 

EtHl  4 

2860 

994 

1.44 

1.20 

1 

39 

EtHl  5 

1264 

593 

0.84 

0.78 

1 

34 

EtHl  6 

2412 

682 

8.70 

4.22 

7 

190 

EtHl  7 

2124 

663 

11.80 

7.74 

5 

192 

EtHl  8 

1586 

567 

23.60 

10.09 

4 

66 

EtHl  9 

2102 

626 

14.70 

10.86 

4 

190 

EtH20 

2204 

646 

5.32 

2.86 

0 

40 

EtH21 

722 

176 

102.89 

68.14 

2 

27 

EtH22 

1848 

607 

17.52 

10.75 

5 

177 

EtH23 

1551 

486 

35.28 

15.10 

2 

89 

EtH24 

701 

165 

171.73 

55.19 

0 

11 
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3.3  f-decalin  +  methyl  radical 

The  third  application  we  considered  was  the  reactivity  of  a  /-decal in  molecule  in  presence  of  a 
methyl  radical.  The  former  is  used  in  surrogate  fuels  and  the  latter  is  a  common  species  in 
combustion  environments  and  plays  an  important  role,  among  other  things,  in  H  abstraction 
reactions.  This  system  can  evolve  in  eleven  different  products  as  shown  in  Figure  3.7,  each 
one  characterized  by  different  multiplicity,  barrier  height  and  shape. 


Figure  3.7:  Possible  reactions  of  the  /-decal in  +  CH3  system.  Six  pathways  involve  the 
unassisted  (1-3)  and  assisted  (la-3a)  C— H  bond  breaking,  four  (4, 5, 7, 8)  the  C— C  cleavage, 
and  one  (6)  the  isomerization  to  c-decalin. 


3.3.1  fe  calculations 

As  before  to  compare  the  rate  computed  with  FERN  we  had  to  first  compute  the  FE  profile 
associate  with  each  reaction.  For  this  purpose  the  reactions  were  separated  in  4  groups:  the 
C— H  bond  breaking  (reactions  1,  2,  and  3),  the  C— C  bond  breaking  (reactions  4,  5,  7,  and  8), 
the  isomerization  (reaction  6)  and  the  hydrogen  abstractions  (reactions  la,  2a,  and  3a).  As 
before,  we  computed  the  FE  at  five  different  temperatures  (among  the  set  of  500,  700,  750,  800, 
1000,  1500  and  2000  K)  and  then  fitted  the  results  with  and  exponential  function.  Simulations 
settings  are  similar  the  the  one  listed  for  ethane  in  section  3.2.1.  A  few  examples  of  the  FE 
profiles  for  bond  breaking  reactions  are  shown  in  Figures  3.8  and  3.9. 

For  the  first  two  groups,  as  shown  in  Figures  3.10  and  3.11  the  differences  between  each 
single  reactions  are  in  general  smaller  than  the  accuracy  of  the  FE  calculations.  Therefore 
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Figure  3.8:  Comparison  of  the  FE  profiles  for  selected  C— C  and  C— H  bond  breaking  at  500 
K;  reactions  are  labeled  according  to  Figure  3.7. 


Figure  3.9:  Comparison  of  the  FE  profiles  for  C— C  bond  breaking  at  1500  K;  reactions  are 
labeled  according  to  Figure  3.7. 


we  used  them  together  to  build  a  general  fitting  function  by  computing  for  each  type  of  bond 
breaking. 

The  hydrogen  abstractions  were  instead  studied  by  simultaneously  biasing  the  specific 
C— H  bond  and  the  distance  between  the  hydrogen  and  the  methyl  radical  carbon.  The  results 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


3.3  f-decalin  +  methyl  radical 


25 


Figure  3.10:  Computed  FE  for  the  C— C  bond  breaking  for  /-decal in  at  different  temperatures. 
Results  are  shown  with  the  95%  confidence  interval.  The  black  line  shows  the  exponential 
fitting. 


Figure  3.11:  Computed  FE  for  the  C— H  bond  breaking  for  /-decal in  at  different  temperatures. 
Results  are  shown  with  the  95%  confidence  interval.  The  black  line  shows  the  exponential 
fitting. 


of  these  simulations  are  reported  in  Figures  3.12,  3.13,  and  3.14. 

The  surprising  result  is  that  there  is  no  communication  between  the  two  basins  correspond¬ 
ing  to  the  stable  states  before  and  after  the  H  abstraction.  Therefore,  the  reactions  la,  2a,  and 
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Figure  3.12:  FE  projections  for  reaction  la  (Figure  3.7)  as  a  function  of  the  distance  of  the 


hydrogen  from  the  C  on  the  decalin  and  the  carbon  on  the  methyl  radical.  Energy  is  reported 
in  kcal/mol  and  isolines  are  drawn  every  25  energy  units. 

3a  are  in  this  FF  effectively  two  step  reactions  in  which  the  presence  of  the  methyl  radical 
does  not  facilitate  the  C— H  bond  breaking.  Alternatively,  if  we  consider  them  as  single  step 
reactions  the  abstraction  must  overcome  a  barrier  of  at  least  50  kcal/mol  higher  than  that 
needed  for  the  unassisted  C— H  bond  breaking  (depending  on  the  assumed  reaction  path).  As 
a  consequence  we  don’t  expect  to  observe  these  reaction  in  a  sample  of  a  few  hundred  FERN 
simulations. 

This  result  is  a  consequence  of  the  choice  of  the  specific  version  of  the  AIREBO  FF,  which 
does  not  change  the  charge  distribution  as  a  function  of  the  distance  of  the  methyl  radical.  A 
more  recent  version  of  this  FF  [24]  or  different  reactive  potentials  may  give  a  more  accurate 
picture  of  this  system’s  reactivity.  Nonetheless,  as  explained  in  the  methodology  section,  the 
validity  of  our  approach  is  not  affected  by  the  choice  of  the  potential  and  as  long  as  the  FERN 
produces  results  consistent  with  the  underlying  FE,  we  can  consider  the  technique  successful. 
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Figure  3.13:  FE  projections  for  reaction  2a  (Figure  3.7)  as  a  function  of  the  distance  of  the 
hydrogen  from  the  C  on  the  decalin  and  the  carbon  on  the  methyl  radical.  Energy  is  reported 
in  kcal/mol  and  isolines  are  drawn  every  25  energy  units. 
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Figure  3.14:  FE  projections  for  reaction  3a  (Figure  3.7)  as  a  function  of  the  distance  of  the 
hydrogen  from  the  C  on  the  decalin  and  the  carbon  on  the  methyl  radical.  Energy  is  reported 
in  kcal/mol  and  isolines  are  drawn  every  25  energy  units. 
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3.3.2  Results 

For  the  FERN  simulations  of  these  systems,  we  followed  the  same  general  protocol  used  for 
ethane,  with  the  addition  of  a  soft  wall  placed  on  the  distance  between  the  center  of  mass  of 
the  methyl  and  the  decalin  molecules  at  0.8  nm.  This  constraint  was  added  to  increase  the 
number  of  collisions  between  the  radical  and  the  decalin  so  that  the  collision  frequency  would 
not  be  a  bottleneck  in  the  assisted  hydrogen  abstraction  rates.  We  tested  different  conditions 
as  listed  in  Table  3.7 

Table  3.7:  List  of  the  different  systems  used  to  test  the  /-decalin  reactivity.  Bias  shape 
parameters  (//#,  Oh)  as  well  as  FE  barriers  (AG*)  are  expressed  in  kcal/mol,  while  bias 
addition  frequency  Sh  is  expressed  in  ps.  In  all  cases  META  algorithm  was  used  (bf  =  °°). 
Reactions  are  marked  according  to  the  labels  defined  in  Figure  3.7. 


Label 

SH 

oH 

hH 

Runs 

1 

2 

3 

AG* 

4  5 

6 

7 

8 

Setl 

0.05 

0.1 

0.1 

181 

104 

104 

104 

108 

108 

115 

108 

108 

Set2 

0.01 

0.1 

0.1 

183 

104 

104 

104 

108 

108 

115 

108 

108 

Set3 

0.05 

0.2 

0.2 

187 

104 

104 

104 

108 

108 

115 

108 

108 

The  results  are  reported  in  Table  3.8  and  the  relative  importance  of  the  different  reactions 
rate  is  compared  to  the  theoretical  value  in  Figure  3.15. 


Set  1  Set  2 

Figure  3.15:  Comparison  between  the  observed  and  theoretical  probability  of  reactions  for  the 
/-decalin  +  CH3  system  for  two  different  set  of  conditions.  Reactions  are  labeled  according  to 
Figure  3.7.  In  both  cases  the  hydrogen  abstractions  (reactions  la  to  3a)  are  not  reported  due 
their  negligible  probability.  Vertical  bars  represent  the  95%  confidence  intervals. 
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Table  3.8:  List  of  the  results  for  the  /-decal in  reactivity  tests.  Temperature  average  and  standard 
deviation  ((T),  cy)  are  expressed  in  K;  time  average  and  standard  deviation  ((/),  at)  are  in  ps. 
Reactions  are  marked  according  to  the  labels  defined  in  Figure  3.7. 


Label 

IT) 

Oj 

(0 

ot 

1 

Occurrences 
2  3  4 

5 

6 

7 

8 

Setl 

2424 

435 

11.52 

7.20 

62 

59 

10 

7 

20 

1 

18 

4 

Set2 

2246 

931 

7.40 

7.08 

64 

61 

14 

6 

14 

3 

11 

10 

Set3 

1928 

1418 

1.27 

1.72 

85 

65 

23 

2 

6 

0 

5 

1 

The  results  show  an  excellent  agreement  for  all  the  pathways  and,  as  before,  we  observe  a 
tremendous  speedup,  with  average  simulation  times  in  a  range  of  2  to  30  ps  depending  on  the 
biasing  parameters.  As  expected,  in  the  FERN  simulations  the  hydrogen  abstraction  reactions 
are  not  observed.  Since  these  are  indicated  in  the  literature  as  the  dominating  pathways  in  such 
conditions  [25],  is  clear  that  the  choice  of  FF  plays  a  crucial  role  in  the  determinations  of  the 
pathways.  In  this  regard,  the  substantial  speed  up  of  FERN  with  respect  to  other  methods  allows 
using  more  accurate  techniques  than  classic  MD  even  with  todays  computational  resources. 
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A  variety  of  natural  phenomena  comprises  a  huge  number  of  competing  reactions  and  short¬ 
lived  intermediates.  Any  study  of  such  processes  requires  the  discovery  and  accurate  modeling 
of  their  underlying  chemical  reaction  network.  However,  this  task  is  challenging  due  to  the 
complexity  in  exploring  all  the  possible  pathways  and  the  high  computational  cost  in  accurately 
modeling  a  large  number  of  reactions.  Fortunately,  very  often  these  processes  are  dominated 
by  only  a  limited  subset  of  the  network’s  reaction  pathways. 

In  this  work,  we  propose  a  novel  method  with  limited  computational  requirements  that  is 
able  to  identify  and  select  the  key  pathways  of  complex  reaction  networks,  so  that  high-level 
ab-initio  calculations  can  be  more  efficiently  targeted  at  these  critical  reactions.  The  method 
estimates  the  relative  importance  of  the  reaction  pathways  for  given  reactants  by  analyzing  the 
accelerated  evolution  of  hundreds  of  replicas  of  the  system  and  detecting  products  formation. 
Within  the  quite  general  validity  of  Equation  1.1  our  method  identifies  the  subset  of  the 
most  likely  reaction  pathways  with  a  minimal  computational  effort  and  without  assumptions 
on  the  reactions  or  transition  states  or  the  need  to  define  different  collective  variables  for 
each  reaction.  Importantly,  the  method  is  efficiently  iterative,  as  it  can  be  straightforwardly 
applied  for  the  most  frequently  observed  products,  therefore  providing  an  effective  algorithm 
to  identify  the  key  reactions  of  extended  chemical  networks.  We  verified  the  validity  of  our 
approach  on  three  different  systems,  including  the  reactivity  of  /-decal in  with  a  methyl  radical, 
and  in  all  cases  the  expected  behavior  was  recovered  within  statistical  error.  These  tests,  in 
particular  the  decalin  reactivity,  with  eleven  different  pathways  that  include  both  unimolecular 
and  bimolecular  reactions,  show  the  full  potential  of  our  approach:  the  FERN  method  is  able 
to  tremendously  speed  up  the  reactivity  of  gas  phase  reactions  without  requiring  any  previous 
knowledge  of  the  system  reactivity.  The  reproducibility  of  the  results  independently  from 
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the  bias  and  simulation  parameters,  as  well  as  the  generality  of  the  potential  energy/SPRlNT 
combination,  makes  this  method  suitable  for  its  effective  iterative  application  to  complex 
reaction  networks.  Moreover,  the  options  of  varying  the  number  of  runs  for  given  reactants 
and  blocking  specific  reactions  allow  a  very  efficient  exploration  even  of  simple  or  partially 
known  reaction  networks. 

Finally,  this  acceleration-detection  approach  can  be  extended  to  different  classes  of  systems, 
with  minimal  adjustments.  For  example,  while  the  potential  energy  and  SPRINT  combination  is 
a  general  streamlined  choice  for  systems  in  the  gas  phase,  for  other  type  of  reactive  networks, 
like  chemical  reactions  in  solutions,  the  relevant  relaxation  times,  e.g.,  water  reorientation, 
should  be  taken  into  account  so  that  the  behavior  of  the  accelerated  system  is  not  biased  by  a 
specific  initial  configuration. 
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